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1 Introduction 


Advection-diffusion-reaction (ADR) equations are widely used for systems in chem¬ 
ical reaction problems. In the linear convection or advection dominated case, stabi¬ 
lized continuous finite elements and discontinuous Galerkin (DG) methods are ca¬ 
pable of handling the nonphysical oscillations. On the other hand, in the non-linear 
stationary case, the non-linear reaction term produces sharp layers in addition to 
the spurious oscillations due to the convection. Thus, an accurate and efficient nu¬ 
merical resolution of such layers is a challenge as the exact location of the layers 
are not known a priori. In the non-stationary case, the resolution of such layers 
is more critical since the nature of the sharp layers may vary as time progresses. 
Recently, several stabilization and shock/discontinuity capturing techniques were 
developed for linear and non-linear steady-state problems H). In contrast to the 
stabilized continuous Galerkin finite element methods, DG methods produce sta¬ 
ble discretizations without the need for stabilization parameters. The DG combines 
the best properties of the finite volume and continuous finite elements methods. Fi¬ 
nite volume methods can only use lower degree polynomials, and continuous finite 
elements methods require higher regularity due to the continuity requirements. The 
DG method is in particular suitable for non-matching grids and hp (space and or¬ 
der) adaptivity, detecting sharp layers and singularities. They are easily adapted 
locally for nonconforming finite elements requiring less regularity. Higher order 
DG approximation can be easily used by hierarchical bases O, and the conver¬ 
gence rates are limited by the consistency error which makes the DG suitable for 
complex fluid flows. The DG methods are robust with respect to the variation of 
the physical parameters like diffusion constant and permeability. The stability of 
the DG approximation retained by adjusting the penalty parameter to penalize the 
jumps at the interface of the elements. Various choices of the penalty parame¬ 
ter is suggested in the literature l[8l[T0l|TTl|2Tl. A unified analysis of the interior 
penalty DG methods for elliptic PDFs are given in lUJ. Other advantages of the 
DG methods are conservation of mass and fluxes and parallelization. Moreover, 
DG methods are better suited for adaptive strategies which are based on (residual- 
based) a posteriori error estimation. The disadvantages are the resulting larger and 
denser matrices and ill-conditioning with increasing degree of the discontinuous 
polynomials. 

The main tool in this paper is the adaptivity applied to the DG discretized (in 
space) systems using (residual-based) a posteriori error estimates. The aim of the a 
posteriori estimates is to derive bounds between the known numerical solution and 
the unknown exact solution. For elliptic diffusion-convection-reaction equations, 
there are various studies applying adaptivity using a posteriori error estimates in 
the literature l[9l[l2l[T3l[T71, which are well-understood. In this paper, we derive 
a posteriori error estimates for semi-linear ADR problems using the well-know el¬ 
liptic a posteriori estimates, by which in contrast to the standard energy techniques 
we do not need to try to adapt the estimates case by case in order to compare the 
exact solution with numerical solution directly. For this reason, we use the elliptic 
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reconstruction technique in llT4ll which allows us to utilize available a posteriori 
estimates derived for elliptic equations to control the main part of the spatial er¬ 
ror. The idea of the elliptic reconstruction technique is to construct an auxiliary 
solution whose difference to the numerical solution can be estimated by a known 
(elliptic) a posteriori estimate, and the constructed auxiliary solution satisfies a 
variant of the given problem with a right hand side which can be controlled in an 
optimal way. By this way, we are able to obtain results being optimal order in both 
and L°°(L^)-type norms, while the results obtained by the standard energy 
methods are only optimal order in L^(/f^)-type norms, but sub-optimal order in 
L°°(L^)-type norms. In Cl, a posteriori error estimates in the U°{L^) -hL^(//^)-type 
norm are derived for linear parabolic diffusion-convection-reaction equations using 
backward Euler in time and discontinuous Galerkin in space utilizing the elliptic 
reconstruction technique. In this paper, we extend the study in Q in a similar way 
by deriving and implementing a posteriori error estimates in the L°°(L^) +L^{H^)- 
type norm using backward Euler in time and discontinuous Galerkin (symmetric 
interior penalty (SIPG)) in space for the convection dominated parabolic problems 
with non-linear reaction mechanisms. To derive the a posteriori error estimates, we 
use the modification of the robust (in Peclet number) a posteriori error estimator in 
uni for linear steady-state diffusion-convection-reaction equations to the steady- 
state diffusion-convection equations with non-linear reaction mechanisms utilizing 
the elliptic reconstruction technique [1141 . 

Application of the adaptive discontinuous Galerkin methods and a posteriori 
error estimates to the problems in geoscience are reviewed recently in O. Most 
of the applications of DG methods in geoscience concern reactive transport with 
advection O [131 HSl and strong permeability contrasts such as layered reservoirs 
ll^ or vanishing and varying diffusivity posing challenges in computations lITSl . 
The permeability in heterogeneous porous and fractured media varies over orders 
of magnitude in space, which results in highly variable flow field, where the local 
transport is dominated by advection or diffusion lIT^ . Accurate and efficient nu¬ 
merical solution of the ADR equations to predict the macroscopic mixing, anoma¬ 
lous transport of the solutes and contaminants for a wide range of parameters like 
permeability and Peclet numbers, different flow velocities and reaction rates and 
reaction rates are challenging problems lIT^ . In order to resolve the complex flow 
patterns accurately, higher order time stepping methods like exponential time step¬ 
ping methods are used lIT^ . We show here using time-space adaptive first order 
backward Euler and DG in space, the same results can be obtained. 

The rest of this paper is organized as follows. In the next section, we introduce 
the function spaces and related norms which are used in our analysis, and the model 
problem with the assumptions to have unique solution. In Section we give the 
symmetric discontinuous interior penalty Galerkin discretized semi-discrete sys¬ 
tem, and fully-discrete system using Backward Euler in time. The derivation of a 
posteriori error estimator is given in Section]^ first for steady-state problems (with 
proofs), then for parabolic problems utilizing the elliptic reconstruction technique. 
In Section we state the adaptive algorithm procedure using the a posteriori er- 
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ror estimators derived in Section The solution procedure of the fully-discrete 
system by Newton method and the structures of the arising matrix and vectors are 
discussed in Section]^ After demonstrating the performance of the algorithm by 
some numerical studies in Section]^ the paper ends with conclusions. 


2 Model problem 


Let C be a bounded, open and convex domain with boundary dQ.. For a 
Banach space X, define the spaces L^(0, T\X) 


xdt 


\M\lp{0,T-,X) — 11^(011 

|v||L-( 0 ,r;X) = esssup||v(0||x < ^ , 

0<t<T 


^ , for 1 < P < +00 
for p = +00 


Also define the space 

(0, T-X) = {v € l2(0, T;X)\ v, G L^{0, T'X)}. 

We denote by C(0, T ;X) and (0, T ;X), the spaces of continuous and Lipschitz- 
continuous functions v : [0, T] ^X, respectively, equipped with the norms 


l|v||c( 0 ,r;X) = 11^(0 11^ < °° 

^llc0.‘(0,r;X) ='^ax{lkllc(0,r;X)) lk/||c(0,r;X)} < °° 


We consider the system of semi-linear diffusion-convection-reaction equations 

^ - V-(8/Vw/)+A--V w/ + r/(w) =/;• , /=1,2,...,/ (3) 


in X (0,r] for the vector of unknowns u = ..., with appropriate 

boundary and initial conditions. We assume that the source functions fi G C(0, T;L^(Q)), 
and the velocity fields j8/ G C (O, T; either given or computed. For flow 

in heterogeneous media in Section |7.4[ the symmetric dispersion tensors 8/ are 
taken of the form 


D} 0 
0 Df 


with 0 < D} ,D? 1. Moreover, we assume that the non-linear reaction terms 

are bounded, locally Lipschitz continuous and monotone, i.e. satisfies for any 
> 0, ^,^1,^2 ^ the following conditions ll^ 


k;(^)|<G, G>0 

|k/(5l)-n-(52)||L2(n) <^lkl -^2||L2(a), ^>0 

neC^(R+), n(o) = o, r'{s)>o. 
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(4a) 

(4b) 

(4c) 





We also assume that there are /c, > 0 satisfying for / = 1,2,...,/ 

^ ^5 II ft||c( 0 ,r;L°°(n)) ^ ^*^5 (5) 

The first inequality in Q ensures the well-posedness of the problem, and the letter 
is used in the a posteriori error analysis. 

The weak formulation of the system ^ reads as: for any v E find 

/= 1,2,...,/, such that for alU E (0,/] 


C du 

/ -^vdx-\-a(t;ui^v) -\-bi(t;u^v) =//(v) 

Ja 

(6) 

a{t;ui,v)= / (8/VwrVv-hA*Vw,-v)/x, 

JQ 

(7a) 

II 

(7b) 

II 

E 

(7c) 


which have unique solutions in the space C{Q^T\L?{Q)) under the given 

regularity assumptions and the conditions (|^-(|4cl). 

In the sequel, for simplicity, we just consider a single equation of the system Q 
(J=l) without any subscript to construct the discontinuous Galerkin discretization 
in the next section and thereafter, and we also continue with a homogeneous dis¬ 
persion tensor leading to a simple diffusivity constant 0 < 8 <C 1. We, furthermore, 
take into account the homogeneous Dirichlet boundary conditions to simplify the 
notations. It can be proceeded with heterogeneous dispersion tensor and other type 
of boundary conditions in a standard way. 

3 Discontinuous Galerkin discretization 

For space discretization of a single equation of Q, we use the symmetric discon¬ 
tinuous interior penalty Galerkin (SIPG) method dlTbl with the up winding for the 
convection part Ena. 

Let {(§/,} be a family of shape regular meshes with the elements (triangles) 
Kj G satisfying Q. = and KjilKj = 0 for Ki, Kj G Let us denote by 
and the set of interior and Dirichlet boundary edges, respectively, so that 
r = U forms the skeleton of the mesh. For any K G let be the set 

of all polynomials of degree at most k on K. Then, set the finite dimensional space 

Vhi^h) = {ve L\a) : G P*(/r), V/f G ^h} t 

The functions in in contrast to the standard (continuous) finite elements, 

are discontinuous along the inter-element boundaries causing that along an interior 
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edge, there are two different traces from the adjacent elements sharing that edge. 
In the light of this fact, let us first introduce some notations before giving the SIPG 
formulation. Let Ki, Kj G (i < j) be two adjacent elements sharing an interior 
edge e = KiDKj cTq (see FiglT]). We denote the trace of a scalar function v from 
inside Kt by vt and from inside Kj by vj, and then we set the jump and average 
values of v on the edge e 

[v] = ViUe-Vjtie, {v} = i(v, + V^). 

Here denotes the unit normal to the edge e oriented from Kt to Kj. Similarly, we 
set the jump and average values of a vector valued function ^ on e 

[q\ = qi ■ He - qj -tie, {^} = ^ + ^j) > 

Observe that [v] is a vector for a scalar function v, while, [^] is scalar for a vector 
valued function q. On the other hand, along any boundary edge e — Kt^] dQ., we 
set 

[v] = {v} = Vi, [q\ = qi ■ n, {q} = qi 

where n is the unit outward normal to the boundary at e. We define the inflow and 
outflow boundary parts at a time t by 

r- = {xeda\p{x,t) ■ n(x) < 0}, r+ = da\r- 

At a time t, the inflow and outflow parts of an element K are defined by 

dK- = {JC e dK\ P{x,t)-nKix) < 0} , dK;^ = dK\dK- 

where H^(x) denotes the outward unit normal to the boundary of the element K at 

X. 



Figure 1: Two adjacent elements sharing an edge (left); an element near to domain 
boundary (right) 

Under the given definitions, the semi-discrete problem reads as: for ^ = 0, set 
Uh{0) G Vh{^h) as the projection (orthogonal -projection) of uq onto Vh{^h)\ for 
each t G (0, T], for all G V/z(<^/z), find G (0, T\Vh{^h)) such that 

C dll 

J^-^Vhdx + ah{f,Uh,Vh) + Kh{uh,Vh) + bh{f,Uh,Vh) = lh{vh), ( 8 ) 
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(9a) 


ah{f,Uh,Vh) = Y, / £'^Uh-yvhdx+ Y / ^■'^UhVhdx 

Ke^h Ke^h 

+ Y [ l-MuT-UHHds-Y ! P 


+ 1 ? / m -[* 

Je 

Kh{uh,Vh) = -Y i{^'^^h}-[uh] + {£'^Uh} ■ [Vh])ds 

bh{f,Uh,Vh) = Y / r{uh)vhdx, 

Ke^kJ’^ 

k{Vh)= Y / fh'^hdx 

k^e,. Jk 


■ UKUhVfids 


(9b) 

(9c) 

(9d) 


where denotes the value on an edge from outside of an element K. The pa¬ 
rameter a G Mq is called the penalty parameter which should be sufficiently large; 
independent of the mesh size h and the diffusion coefficient e ifT^ [Sec. 2.7.1], We 
choose the penalty parameter a for the SIPG method depending on the polynomial 
degree has <y — 3k{k+ 1) on interior edges and c = 6k{k+ 1) on boundary edges. 

Upon integration by parts on the convective term, bilinear form ( |9^ will be 

ah{f,Uh,Vh) = Y / £^Uh-yvhdx- Y / {^Uh-^Vh + y •^UhVh)dx (10a) 

K€^h Ke^h 

+ Y p ■nKUh{vh-vl'“)ds+ Y / P-HKUhVhds 

+ Y^ lH]-[Mds 

.tf Je 


Note that the bilinear form ah{t', w, v) is well-defined for the functions w, v G /^q (f2), 
and equal to 


ah{t\u^v) — / (8Vw• Vv +j8 • Vmv)Jx 

Ja 

Thus, the weak formulation can be rewritten for any ^ G (0, T] as 
C du 

/ —vdx + ah{t\u^v)+b{t\u,v) = /(v) , Vv G 

Ja 


( 11 ) 


For the fully discrete system, consider a subdivision of [0,7] into n time in¬ 
tervals 4 = of length Tk, k— 1,2,... ,az. Set = 0 and for — 

^1 + ^2 H-Denote by an initial mesh and by the mesh associated to the 

time step for /: > 0, which is obtained from by locally refining/coarsening. 
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Moreover, we assign the finite element space Vj^ = Vh{^^) to each mesh Then, 
the fully discrete problem, backward Euler in time, reads as: for ^ = 0, set G 
as the projection (orthogonal -projection) of uq onto for ^ = 1,2,... ,n, find 
ul G Vj^ such that for all G Vj^ 




4-^h 


k-\ 


% 


-vldx + ah{t'^\ul,\i)+Kh{ul,\i)+bh{t'^-,ul,vl)= fv^dx ( 12 ) 

Ja 


4 A posteriori error analysis 


The problems concerned in this paper are convection/reaction dominated non¬ 
stationary models which often produce internal/boundary layers where the solu¬ 
tions have steep gradients. Thus, an efficient solution of such equations are needed. 
The most popular technique to solve the convection/reaction dominated problems 
is the adaptivity which requires an estimator to locate the regions where the error is 
too large. In this work, we construct the residual-based, robust (in Peclet number) 
a posteriori error estimators for non-stationary convection domianted diffusion- 
convection equations with non-linear reaction mechanisms. To construct the a pos¬ 
teriori error estimators, we extend the a posteriori error estimators for linear non¬ 
stationary problems constructed in [|71 which uses the a posteriori error estimators 
for linear stationary models constructed in ifTTll utilizing the elliptic reconstruction 
technique lfT4l to make connection between the stationary and non-stationary error. 
As in lO, first we construct and prove the a posteriori error bounds for non-linear 
stationary models in Section |4T| Then, in Section|4^ we give the a posteriori error 
bounds for the semi-discrete system of the non-stationary problems with non-linear 
reaction term. Finally, in Section [43| we give the a posteriori error bounds for the 
fully-discrete system of the non-stationary problems with non-linear reaction term. 
The main contribution to the error analysis lies in the construction of the error 


bounds for non-linear stationary models in Section [4T| The remaining utilize the 
constructions in ||7]|. 


4.1 A posteriori error bounds for non-linear stationary system 

We consider first the stationary problem, i.e. for a given t G (0, T], let G Hq (Q.) 
be the unique solution to the weak problem 

a{t\u\v) +b{t\u\v) = /(v) , Vv G (13) 

and let u\^Vh{^h) be the unique solution to the discrete problem 

= ^ VvGy/,((^/,) (14) 

In order to measure the spatial error, we use the energy norm 

p 

llkllP= T (ll^^^lli2(;f)+'^ll^llL2(^))+T'^llMlli2(e)) 

K€^h eer 
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and the semi-norm 


where 


eer ^ 

, , Tq w • Vwdx 

|w|* = sup —, 

w€Ho‘(n)\{o} infill 

For each element K e ^h,'we define the local error indicators as 


(15) 


where 

,2 




^Rk ~ PK\\fh +fiK)\\h(^lC), 


(16) 


nL = 


Pjk = 




' e€dK\Y^ 


1 ^ ^ EG 

^ ^ KJTe 


' e€dK\Y^ 


, K 

+ Khg H- 


■SI Il 2 


eedKCY^ 


Z.2(e)+ L \ 


..s l|2 


(17) 


with the weights Pk and pg, on an element K, given by 

P/f = min{/z^£"2,K-“2}, Pe = min{fi^e"5,K:"2}, 

for K" 7 ^ 0. When fc = 0, we take pK = hK£~^ and Pe = /?<,£“5. Then, our a 
posteriori error indicator is given by 

^ = f E 

\Ke^k ) 

We also introduce the data approximation terms, 

®i ^ P1(II/~ A|Il 2(/^) + ll(^ ~Ph) • ^K\\h{K))‘ 

and the data approximation error 

• (19) 

\Ke^h J 

Then, we have the a posteriori error bounds 


— ^h\\DG ^ t] +0 

^^II«'-«aIIdg + © 


(reliability), 
(efficiency). 


( 20 ) 

( 21 ) 


with 
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4.1.1 Proof of a posteriori error bounds 

The proof for the bounds of spatial error estimates are analogous to the ones in fTll 
for the linear problems. Therefore, only the proofs for non-linear reaction term are 
stated explicitly . In the following, we use the symbols < and > to denote the 
bounds that are valid up to positive constants independent of the local mesh size h, 
the diffusion coefficient £ and the penalty parameter a. 

The spatial error \\u^ — is not well-defined, since G and G 

Vh{^h) ^ ^0 (^)- Therefore, we split the stationary SIPG solution as 

with G hI {Q)r\Vh{^h) being the conforming part of the solution and G Vh{^h) 
is the remainder term. In this way, we have G +Vh{^h)^ and from the 

triangular inequality 


\W — '^AiWdg < \W — '^AiWdg^ II^^I|£>g 

Now, both the terms on the right hand side are well-defined norms, and our aim is 
to find bounds for them. Next, we introduce the following auxiliary forms: 


Ke^h 


D/^(w,v) = ^ / TeVw • Vv — V • j8wv^ Jx 


• riKUvds 


Oh{u^v) —— ^ f Pu’Vvdx-\- ^ f P 

-h ^ P -nKU^v — v^^^)ds 

Ke^JdK+\dQ 

Jh{u,v) = Y, ^ 

eGFoUr Je 

We note that, for a specific f G (0, T], the SIPG bilinear form (|10a) satisfies 


(22a) 


(22b) 

(22c) 


ah{f,u,v) = Dh{u,v) + Oh{u,v) + Jh{u,v) 

and is well-defined on Hq (Q) -|-V/,(<§/,). Using the first identity in we can easily 
have for any u E Hq (Q) 

afi(t;u,v) > |||m||P 

Moreover, the auxiliary form s are continuous ifTTI ILemma 4.2]: 


\Dh{u,v)\ < |||u||| |||v 

\Oh{u,v)\ < |^m|* |||v| 

l•4(M,v)| < |||n||| |||v 


u,v(^Hl{a) + Vh{^h), (23) 

M G //J (a) + 14 (<§;,), V G Hi (Q), (24) 
M,vG//J(f2) +14 (<§/,), (25) 


10 





and for u e Vh{^h), v G Vh{^h) nf/d(Q) ifTTl [Lemma 4.3] 


\Kh{u,v)\ < £ 

V^GroUi 




1/2 


(26) 


VeGFoUrO « 

We also have for the non-linear form bh{t;u,v), for a specific time t, using the 
boundedness assumption ( |4a| ), 

\h{f,u,v)\ < |||v|||, u,v EHQ{Q.) + Vh{^h)- (27) 

Now, we give some auxiliary results and conditions which are used in the proofs. 

• Inf-sup condition: lUTl [Lemma 4.4] For a nonzero u G Hq{£1), for a con¬ 
stant C > 0, we have 


|m||| -f \Pu\ 


< 


sup 


ah{t;u,v) 


(28) 


• Approximation operator: Let = Vh{^h) (^) the conforming sub¬ 
space of F/j((§/j). For any m G F/j (<§/,), there exists an approximation operator 
Ah : Vhi^h) ^ satisfying 

£ < £ f he\[u]f-ds, (29) 

Ke^ eeFoUTD-'e 

Y,\\'^iu-Ahu)\\l2(^K)< £ [ ^\[u]\^ds. (30) 

Ke^ eGFoUrO 7e 


• Interpolation operator: For any u G Hq{Q.), there exists an interpolation 
operator 

4 : ^ {wGC(Q) : w|/f G Pi( 7 :),V 7 : G = 0 on F} 


that satisfies 


\\\hu\ 


< 


1/2 


£ Pk 11^ ) 

\Ke^ J 

£ e^'^Pe^h-huWli^K)] 

VeGFoUrO / 


< 


1/2 


< 


(31) 

(32) 

(33) 


Now, consider the splitting of the stationary solution ul = ul + u'j^ as = 
AhU^ G Hq{Q.) nVh{^h) with Ah is the approximation operator and m]" = 
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• Bound for remainder term: IITtI rLemma 4.7] For the remainder term, we 
have the bound 

hhWDG^ri (34) 

where 7] is our a posteriori error indicator ( [T^ . 


Lemma: For a given 7 G (0, T] and for any v gHq (tl), we have 

/ f{v-Ihv)dx-ah{t;ul,v-Ihv)-bh{t;ul,v-Ihv) < (t? + ©)||h 
Ja 

where 4 is the interpolation operator. 

Proof: Let 

^ = / f{v-Ihv)dx-ah{t\ul,v-Ihv)-bh{t\ul,v-Ihv). 

Integration by parts yields 

^ = L / if + ^^^h-P-^K-riul))iv-Ihv)dx 
Jk 

sVu^ ■ nK{v — Ihv)ds 


Ke^h 

I 


+ Y. i ^■nk{ul-Uif"^){v-Ihv)ds 

= T 1 +T 2 + T 3 . 

Adding and substracting the data approximation terms into the term Ti 

T\ ^ Y [ {fh + £^ul-'^h-^ul-r{ul)){v-Ihv)dx 


Y [ (if-fh)-{P-Ph)-^lll){v-Ihv)dx. 


(35) 


Using the Cauchy-Schwarz inequality and interpolation operator identity ( [32| ) 

1/2 / 1/2 

KKe^h 


Ke^H J KKe^h ) 

Y PK'^W^-^hvWll^K) 

XKe^h / 


< 


-Lei 

L (^1;.+®: 

KKe^h 


1/2 
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For the terms Tz and we have IITtI [Lemma 4.8] 


T2<(Lnl) 


1/2 


1/2 


7^3 < I 4 


KKe^h 


which finishes the proof. 


Bound to the conforming part of the error: For a given ^ G (0, T], there holds 

\W-ul\\DG<ri + e. (36) 

Proof: Since — Hq (H), we have |m'' — |c = |j3 («'' — m^) |*. Then, from the 


inf-sup condition (28) 


\u^ — M^IIdG = — M^lll + !«■' — M^lc ^ sup 

veHi{a)\{ 0 } 


ah{t;u^ - ul,v) 


So, we need to bound the term ah{f,u^ — m^,v). Using that — u'j^ ^ Hq{Q.), we 
have 

ah{t;u^-ul,v) = ah{t;u\v)-ah{t;ul,v) 

= / fvdx-bh{f,ufv)-ah{f,ul,v) 

Ja 

= / fvdx-bhit;u\v)-Dh{ul,v)-Jhiul,v)-Ohiul,v) 

JQ. 

= / fvdx-bh{f,ul,v)+bh{f,ul,v)-bh{f,ufv) 

Ja 

-ah{t;ul,v) + Dh{ul,v) + Jh{ul,v) + Oh{ul,v). 

We also have from the SIPG scheme 

/ fhvdx = ah{t\ulfhv)+Kh{ul,Ihv)+bh{t-,ul,Ihv) 

Ja 


Hence, we obtain 


a{t;-ul,v) = Ti + T 2 + T 3 + T 4 


Tx = 


/ f{v-Ihv)dx-ah{r,uf,v-Ihv)-bh{t\uf,v-Ihv) 
Ja 

T2 = Dh{ul,v)+Jh{ul,v) + Oh{ul,v) 

T3 = Kh{ul,Ihv) 

T 4 = bh{f,ul,v)-bh{f,ufv) 


13 




From the inequality ( [35] ), we have 

Ti<{n + ®)\\\v\\\ 

The continuity results ( |23||25] ) and the bound to remainder term ( [3^ yields 

72<(||KIII + IMI*)|||v|||<t^ 

Moreover, using the identities ( [26| ) and ( [3T] ), we get 


Ve^ ) \Ke^ 


1/2 


Finally, using Cauchy-Shwarz inequality and the boundedness property ( |27) ), we 
get 

T4 = bh{t\u\^v)—bh{t\u^^v) — / r{ul)vdx— r{u^)vdx 

Jo, Jo, 

^ ^l||^llL2(n) ~ ^ 2 \\ v \\ i 2 (^ q ^ 

This finishes the proof. 

Proof to the reliability: Combining the bounds ( [34| ) and ( |36| ) to the remainder and 
the conforming parts of the error, respectively, we obtain 

IW^ — ^hWoG < 11^*^ ~ W^||z)G+||w^||z)G 

< 77+0 + 77 

< T7+© 

Proof to the efficiency: The proof of the efficiency is similar to Theorem 3.3 in 
113. We only use the boundedness property ( |4a| ) of the non-linear reaction term to 
bound the terms occurring in the procedure in fTTi . 


4.2 A posteriori error bounds for the semi-discrete system 

In order to measure the error for the semi-discrete problem, we use the L°°(L^) + 
- type norm 


(0,7+2 (n)) 


+ 


^dt 


We make use the elliptic reconstruction technique in llT4l : the elliptic reconstruc¬ 
tion w gHq (Q.) is defined as the unique solution of the problem 


a{t',w,v)-\-bh{t',w,v) = 




^]vdx, \/veHQ{Q.). 


(37) 
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The SIPG discretization, on the other hand, of the above system reads as: for each 
t G (0, T], find w* G C®’* (0, T;Vh{^h)) such that for all v/, G Vh{^h) 


ah{t;wh,Vh) + Kh{wh,Vh) +bh{f,Wh,Vh) = 
which implies that Wh = Uh. Hence, the error bound to the term \\w — Uh\\DG can be 


found using the a posteriori error bound (20) for non-linear stationary problem. 

We give the a posteriori error bounds for the semi-discrete system developed 
as in Q. For the error e{t) = u{t) — Uh{t) of the semi-discrete problem, we set the 
decomposition e{t) = ji{t) + v(/) with ii{t) — u{t) — w{t) and V = w{t) — Uh{t). 
Then, for any t G (0, T], using and ( [T^ , we obtain that 


^vdx-\-ah{t‘ld,v)+bh{t‘lJ.,v)=() 


which yields the error bound fTl [Theorem 5.4] 


< 


where the error estimator fj is defined by 


V = Ik 

with 


= Y. pk 

Ke^h 


T ( / T T ^ 


f - + £Am* - P • Vuh - r{uh) 


L^{K) een 


eerour® 


+ E ( + + ) IIK]llL2(e) 


duh 


il = E 

eGFoUrO 

vl = E ^^IINIIl2(, 

eGFoUrO 


L2(e) 


and the weight p 7 = min(£ 2 ^K 2 ). 


4.3 A posteriori error bounds for the fully-discrete system 

For the fully discrete case, we consider the solutions at discrete time instances. For 
this reason, let A* G Vj^ be the unique solution of the stationary system 

ah{t''-,ui,vl) + Kh{ui,vl)+bh{t''-,ui,vl)= [ A’^v\dx , Vv^ G 

JQ. 
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Note that for k> 1, — {u\ — j with being the L^-projection 

operator onto the space V]^. Then, the elliptic reconstruction G Hq (Q.) is defined 
as the unique solution of the stationary system 

+ Z7(/^;w^,v) = f A^vdx ^ VvG//d(t2). (38) 

Jo. 

Next, we take, as in ||71, the discrete solution Uh{t) as a piecewise continuous 
function so that on each interval Uh{t) is the linear interpolation of the 

values and given by 


Uh{t) = ^+ 4(^)4 

with the linear Lagrange interpolation basis 4-i and 4 defined on Then, 

for the error e — u — Uh, using (11) and (38), we obtain for all v G /^q (f2) 


J^^vdx + ah{t;e,v) + bh{f,e,v) = 

-ah{f,Uh,v) - bh{f,Uh,v) 


duh 


vdx 


which yields the error bound Q [Theorem 6.5] 

Ikll* <T]| + T?r 

where r\s is the spatial estimator given by lITI 


with 




1 ^ 




k=l 


k=l 


+ + min ] ( £ ) ,p^ £ rjill 


<k^l 


k^l 


£ pI X^ + eMl-p-Vul-riU 

+ ^ II [“ft] Ili2(<;) 

I pI ^ ^ 


lHk) 


+ £ e ^^^Pe||[eVM^]||^ 2 (, 


eSFo 


e) 


(39) 


'Ck 


L^{K) 


I/^.ii[4]iiL 


eeT 

eer 


e) 


“a 


% 


L\e) 


16 




















and riT is the temporal estimator given by f7l 




n rtk \ I n rtk \ n rtu 

L / + £ / riT2,kdt] ,PtY. 

k=l I \A:=1 ^h-\ / k=l ^h-\ 


nl 


dt 


(40) 


with 

- M + 4-1 - M~' whiK) 

= Z ll/-/ + 4-i(A^-A^-') + 4(V-i8^-V-i8)4 + 4(V-i8^--i-V-i^^^^^ 


5 Adaptive algorithm 


The time-space adaptive algorithm for linear convection-diffusion problems in Q 
is modified for semi-linear problems of type ([T]), see Fig. It is based on the 
residual-based a posteriori error estimators given in the previous sections. The al¬ 
gorithm starts with an initial uniform mesh in space and with a given initial time 
step. At each time step, the space and time-step are adaptively arranged according 
to the user defined tolerances ttol for time-step refinement, and stol+ and stol“ for 
spatial mesh, former corresponding to refinement and latter corresponding to the 
coarsening procedures in space. Note that we do not need a temporal tolerance 
corresponding to time-step coarsening, since we start, in our problems, with a uni¬ 
form equispaced distribution of [0,r] having sufficiently large time-steps. Thus, 
it is enough just bisecting the time intervals not satisfying the temporal tolerance 
ttol. Both the refinement and coarsening processes in space are determined by the 
indicator rfsik appearing in the spatial estimator. Since the temporal estimator 
riT (|4^ is not easy to compute, the adaptive refinement of the time-steps are driven 
by the modified temporal-indicator [|71 

rtk rtk 

tin = / + min{pr ,T} 

-'4-1 ’ ^tk-i 


sum of which gives a bound for the temporal estimator T]^. 

Although the adaptive algorithm, Fig.[^ stands for a single equation of the sys¬ 
tem (01, it is not difficult to extend the algorithm to the coupled systems. For this, 
say we have a system of two equations for the unknowns ui and U 2 , the temporal- 
indicators ffl , fjr and the spatial-indicators 7] I ,7]? corresponding to the un- 
knowns ui and U 2 , respectively, are computed. To adapt the time-step size, we 
ask the temporal condition for the both temporal-indicators, i.e. 
f\j^ < ttol. On the other hand, to select the elements to be refined, we take the set 
of elements which is the union of the sets of the elements satisfying > stol+ 
and Vsik^ stol+, and similar procedure to select the elements to be coarsened, but 
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Figure 2: Adaptive algorithm chart 








not including any elements which are selected to be refined. Numerical studies 
demonstrate that the adaptive algorithm is capable of resolving the layers in space 
as the time progresses. 

6 Solution of the fully-discrete system 

In this section, we discuss the solution of the fully-discrete system ( p^ on an 
arbitrary time-step, which is solved for all ^ = 1,2,... ,fi. In order to not be 
confused about the notations, let us consider the system ( [T^ on an arbitrary 
time-step without the superscript for the time-step of the form 


/ ——'^Vhdx + ah{uh,Vh)+Kh{uh,Vh)+bh{uh,Vh) ^ / fvhdx, Vv/, G 14 
Ja ^ Ja 

(41) 

where we set Uh := Wh := Vh := v^, / := T := T^, ah{uh,Vh) := v^), 

bh{uh->^h) •= Vh- The approximate solution and the 

known solution (from the previous time-step) Wh of © have the form 


NelNloc NelNloc 

M/I = Xm ’ ^/i = L L ^‘i^i 

1=1 i=l 1=1 


(42) 


where 0/’s are the basis polynomials spanning the space \4, wj’s are the unknown 
coefficients to be found, wj’s are the known coefficients. Nel denotes the num¬ 
ber of triangles and Nloc is the number of local dimension depending on the 
degree of polynomials k (in 2D, Nloc = (^-h l)(^-h 2)/2). In DG methods, we 
choose the piecewise basis polynomials 0/’s in such a way that each basis func¬ 
tion has only one triangle as a support, i.e. we choose on a specific triangle Ke, 
e G {1,2,... the basis polynomials 0^^ which are identically zero outside 

the triangle Ke, I = 1,2,... ,M(9c. By this construction, the stiffness and mass 
matrices obtained by DG methods has a block structure, each of which related 
to a triangle (there is no overlapping as in continuous FEM case). The product 
dof Nel ^Nloc gives the degree of freedom in DG methods. Inserting the linear 
combinations ( |42| ) of Uh and Wh in ( [4T] ) and choosing the test functions Vh = 0/, 
/ = 1,2,... ^Nloc, / = 1,2,... ^Nel, the discrete residual of the system •ED in ma¬ 
trix vector form is given by 


Res{u) = Mu — Mw -\- t{Su + b{u) —/) = 0 (43) 

where w, w G are the vector of unknown and known coefficients wj’s and wj’s, 
respectively, M G is the mass matrix, S G is the stiffness matrix 

corresponding to the bilinear form dh{uh^Vh) := cLhi^hNh) + Kh{uh^Vh), b G 
is the vector function of u related to the non-linear form bh{uh^Vh) and / G is 
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the vector to the linear form fvhdx. The explicit definitions are given by 


Sn 

Sn ■ 

Sl^Nel 


■ Mil 

Mn ■ 


^21 

S 22 


, M = 

M 21 

M 22 


SNel,l 


SNel,Nel_ 


MNel,l 


MNel,Nel_ 



Ul 


Wi 


■ ^ i ( m ) ■ 


' h ' 

u — 

U2 

, ^ = 

W2 

, b{u) = 

hiu) 

, / = 

h 

dNel_ 

_pNel_ 

bNel{u)_ 

JncL 


where all the blocks have dimension Nloc\ 


4(01,01 ) 4(02, 01 ) 

_4(01,07V/(?c) 


4(0^/<5c’01 ) 

^hi^Nioc-^^Nloc). 


/^0{0/jx /^0l0/^X 

-Iq 01 ^Nloc^^ 


Iq ^Nloc^l^ 

Iq ^NIoc^NIoA^- 


Hi = 

u\ 

u ‘2 

, = 

w\ 

II 

.. 

1 _ 

, fi= 

1 - 

1 _ 


J^'nIoc. 


ymoc. 


_4(w/z, 0^/<5c). 


-In f^Nioc^^^. 


The block structure of DG methods causes to the increase of the condition number 
of the obtained matrices by the degree k of basis polynomials, which is a drawback 
of DG methods comparing to the classical (continuous) finite elements. However, 
taking into account the locality, a valuable property, of DG methods, this drawback 
can be handled by various preconditioners developed for DG discretized schemes 
in the literature. Besides, the condition number of the stiffness matrix S increases 
linearly by the penalty parameter <7, as well. Therefor, the penalty parameter 
should not be chosen too large. On the other hand, it should be selected suffi¬ 
ciently large to ensure the coercivity of the bilinear form ifT^ rSec. 27.1], which is 
needed for the stability of the convergence of the DG method. It ensures that the 
stiffness matrix arising from the DG discretization is symmetric positive definite. 
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In the literature, several choices of the penalty parameter are suggested. In iFTTl . 
computable lower bounds are derived, and in ifTOl . the penalty parameter is chosen 
depending on the diffusion coefficient £. The effect of the penalty parameter on the 
condition number was discussed in detail for the DG discretization of the Poisson 
equation in ||8l and in ll^ for layered reservoirs with strong permeability contrasts. 
In our study, we select the penalty parameter a depending only on the polynomial 
degree as O’ = 3k{k+ 1) on interior edges and a = 6 k{k+ 1) on boundary edges. 
The reason of the coupling of the penalty parameter <7 on boundary edges is suf¬ 
ficient penalization of the solution on the boundary due to the non-homogeneous 
Dirichlet boundary conditions which are imposed weakly in DG methods. 

Next, we solve the system by Newton method. In the sequel, we start 
with an initial guess (most possibly = w, i.e. the known solution from the 
previous time-step) and we solve the system 

= -Resiu^^'>) 

^(^+1) = + 5 = 0,1,... 

where = M + t( 5 + 7|) is the Jacobian matrix of the system at the iteration 
and Ji denotes the Jacobian matrix to the non-linear vector b(u) at the iteration . 

1 Numerical studies 

In this section, we give the numerical studies demonstrating the performance of the 
time-space adaptive algorithm. All the computations are implemented on MATLAB- 
R2014a. In the problems, by the very coarse initial mesh, we mean an initial 
mesh which is formed, for instance on = (0,1)^, by dividing the domain with 
Axi = Ax 2 = 0.5 leading to 8 triangular elements and 48 DoFs for quadratic ele¬ 
ments. As the first example. Example |7.1[ we give a test example with polynomial 
type non-linearity having a non-moving internal layer to figure out the benchmark 
of the algorithm by using different tolerances and diffusion parameters 8: rates of 
error, spatial and temporal estimators, and effectivity indices (proportion of the es¬ 
timator to the error). We expect that the effectivity indices lie in a small band for 
different diffusion parameters meaning that our estimators are robust in the system 
Peclet number. Moreover, to demonstrate the mentioned properties, we use the 
weighted DoFs as in ||7| 


1 ^ 

Weighted DoFs = — ^ 

^ k=i 

where Xk denotes the total number of DoFs on the union mesh U Since the 
first example has a non-moving internal layer, a monotonic increase in the DoFs 
is expected by the time progresses. Conversely, we give problems having moving 
layers by the time progresses in Example |7.2||7.3[ In this case, we expect that 
the refinement and coarsening procedures in space work simultaneously leading to 
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oscillations in time vs DoFs plots. By Example |7.2[ we also test the performance 
of our algorithm for a coupled system. As the final example, Example |7.4[ we 
consider an important real life problem representing a reaction in porous media 
having internal layers at different locations due to the high-permeability rocks. 

7.1 Example with polynomial type non-linearity (benchmark of the 
algorithm) 

The first example is taken from ||4l with a polynomial non-linear term 

-h j8 • Vw — eAw-h r(w) = / in tl = (0,1)^, 

with the convection field j8 (x,y) = (2,3)^, diffusion coefficient £ = 10“^, the non¬ 
linear reaction term r{u) = The source function / and the Dirichlet boundary 
condition are chosen so that the exact solution is given by 

u{x^t) = 16sin(7rt)xi(l —xi)x 2 (l —X 2 ) 

[0.5 + 71-^ arctan(2e-'/2('o.252 - (xi - 0.5)^ - (x2 - 0.5)^))]. 

We start by demonstrating the decrease of the errors by uniform time-space refine¬ 
ment using linear DG elements. In Fig. the expected first order convergence in 
space and time is shown. 


10 ^ 


10 ^ 

10 ° 


Figure 3: Example [7. l[ Decays of estimators and errors for uniform time-space 

For the time-space adaptive solution, we use quadratic DG elements. We in¬ 
vestigate the performance of the spatial estimator by fixing the temporal time-step 
T = 0.005 so that the temporal error dominated by the spatial error. We reduce 
the spatial estimator tolerance stol+ from 10“^ to 10“^. The rates of the error and 
the spatial estimator are similar as illustrated in Fig.|^for £ = 10“^ and £ = 10“^. 
Fig.S shows the spatial effectivity indices and the decrease of the spatial estima¬ 
tors for the diffusion constant £. The effectivity indices do not exceed 7 which are 
acceptable as in Q for linear convection-diffusion problems. 


Error 

-•-Spatial Est. 
-•-Temporal Est. 



Weighted DoFs 
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Weighted DoFs 


Figure 4: Example |7.l[ Error/spatial esti m ator for e = 10 ^ (left) and £ = 10 ^ 
(right) 


Spatial Effectivity Indices 




Figure 5: Example [tT] Spatial effectivity indices (left) and estimators (right) 


To investigate the performance of the temporal estimator, we fix a sufficiently 
fine spatial mesh so that the the spatial error dominated by the temporal error, and 
then we reduce the temporal estimator tolerance ttol in the range 10“^ — 10“^. In 
Fig.|^ the temporal effectivity indices and the decrease of the temporal estimators 
are not affected by 8, and effectivity indices are almost the same within the band 
1 - 2 . 

Finally, we apply the time-space adaptive algorithm with the tolerances ttol = 
10“^, stol+ = 3 X 10“"^ and stol“ = 3 x 10“^. Firstly, we prepare an initial mesh 
starting from a very coarse spatial mesh and a uniform partition of the time interval 
[0,0.5] with the step-size T = 0.25 until the the user defined tolerances ttol and 
stol+ are satisfied. The adaptive mesh at the final time T = 0.5 is shown in Fig.[^ 
In Fig. on the right, the change of the time-steps is shown, whereas the change in 
the DoFs is illustrated in Fig. on the left. Since the layers in the problem do not 
move as the time progresses, the number of DoFs increases monotonically by the 
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Temporal Effectivity Indices 



Rates of Temporal Estimator 



Figure 6: Example |7.l[ Temporal effectivity indices (left) and estimators (right) 


spatial grid refinement. In Fig. it is shown that all the oscillations are damped 
out by adaptive algorithm using less DoFs compared to the uniform one. 


DoFs: 98352 



Figure 7: ExampleAdaptive mesh 


7.2 Coupled example with polynomial type non-linearity 

The next example is a coupled non-linear problem taken from lISl . 


— £AUi-^Pi-VUi-\-UiU2 = //, /= 1,2 

ot 


on = (0,1)^ with the convection fields j8i = (1,0)^ and J82 = (—1,0)^, and the 
diffusion constant 8 = 10“^. The Dirichlet boundary conditions, initial conditions 
and the load functions fi are chosen so that the exact solutions are 


u\(x^t) 


1 

2 


^1 — tanh 


2xi — 0.2^ — 0.8\ 

7^ J 
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DoFs: 196608, t=0.005 


DoFs: 98352 



Figure 8 : ExampleUniform (left) and adaptive (right) solutions at T=0.5 


U 2 {x^t) = - ( 1 + tanh 


2x\ +0.2^ — 0.9 


2 ' . \/^ 

We use again quadratic DG elements. We prepare an initial mesh, Fig. on the 
left, starting with a very coarse spatial mesh and a uniform partition of the time 
interval [ 0 , 1 ] with the step-size T = 0.1 until the user defined tolerances ttol = 
10“^ and stol+ = 10“^ are satisfied. Here, two sharp fronts move towards to each 
other and then mix directly after the time / = 0.1, Fig. The movement of the 
fronts are also visible in Fig.[TT]claiming that refinement/coarsening of the adaptive 
algorithm works well. We see that the sharp fronts in the cross-wind direction 
X 2 — 0.5x1 + 0.75 are almost damped out. Moreover, Fig. [TT][T^ show that both the 
spatial and temporal estimators catch the time where the two sharp fronts mix. 


7.3 Non-linear ADR equation in homogeneous porous media 

We consider the advection-diffusion-reaction (ADR) equation in |[T^ with polyno¬ 
mial type non-linear reaction 

—eAw-h j 8 • Vw-h — 1) = 0 intlx(0,r] 
ot 

on = (0,1)^. We take as in the homogeneous dispersion tensor as 8 = 10“^, 
the velocity field j 8 = (-0.01,-0.01)^ and 7 = 100. The initial and boundary 
conditions are chosen by the exact solution 

u{x^t) = [1 -hexp(a(xi +X 2 — bt) -\-a{b— 1 ))]~^ 

with a = y^77(4^ and b = — 0.02 -h The problem is a transport of a front 
in homogeneous porous media. We simulate the given problem for the final time 
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X 10^ 




Figure 9: Example |7.l[ Evolution of DoFs (left) and time-steps T (right) 



t=i 



Figure 10: Example [7.2[ Cross-section plots in the cross-wind direction at / = 0.1 
(left) and / = 1 (right) 


r = 1, and with quadratic DG elements. We begin by preparing an initial mesh 
starting from a very coarse spatial mesh and a uniform partition of the time interval 
[0,1] with the step-size T = 0.25 until the user defined tolerances ttol = 3 x 10“^ 
and stol+ = 10 “^ are satisfied. In Fig.[^ the adaptive meshes and solution profiles 
are shown at the times t — {0.2,0.6,1} where the movement of the front can be 
seen. The time vs DoFs and time vs time step-size plots in Fig. [^indicate clearly 
the oscillations in DoFs and time-steps due to the movement of the front. 
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Figure 11: Example [t^ Adaptive meshes at ^ = 0, ^ = 0.1 and / = 1 (from left to 
right) 
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Figure 12: ExampleEvolution of DoFs (left) and time-steps T (right) 


7.4 Non-linear ADR in deterministic heterogeneous porous media 

We consider the ADR equation in lIT^ with Monod or Langmuir isotherm type 
non-linear reaction 


(j U V / \ .—7 ^ 

^ — V - (eVu) + P{x) • Vm+ —— 
ot \ u 

u{x^t) 

—£'Vu{x^t) -n 

m(x,0) 


0 in X (0,r] 

1 onr^x[o,r] 

0 on {da\r^)x[o,T] 
0 in Q, 


with = (0,3) X (0,2) and E^ = {0} x [0,2]. The problem represents a reaction 
in porous media, for instance, transport in a highly idealized fracture pattern. Here 
£ stands for the heterogeneous dispersion tensor given by 

_ rio-3 0 

^ “ 0 10-4 
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DoFs: 10452 


DoFs: 14628 


DoFs: 16422 



Figure 13: Example |7.3[ Adaptive meshes (top) and solution profiles (bottom) at 
t = 0.2, t = 0.6 and / = 1 (from left to right) 


The velocity field j8 (x) is determined via the Darcy’s law 

p = -^^Vp 


where p is the fluid pressure, p is the fluid viscosity and k{x) is the permeability 
of the porous medium. Using the mass conservation property V • j8(x) = 0 under 
the assumption that rock and fluids are incompressible, the velocity field j8 (x) is 
computed by solving 



P 

P 

—k{x)Vp-n 


0 in 

1 on {0}x [0,2] 

0 on {3}x [0,2] 

0 on (0,3) X {0,2} 


We simulate the given problem for the final time T — I using linear DG ele¬ 
ments. We take the fluid viscosity /i = 0.1, and the permeability field as in ifT^ 
with three parallel streaks the permeability of which are 100 times greater than the 
permeability of the surrounding domain, see Fig. on the left, by which the flow 
is canalized from the lower-permeability rocks into the high-permeability ones. 
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Figure 14: Example [?. 3 [ Evolution of DoFs (left) and time-steps T (right) 


Fig. on the right. For the adaptive procedure, we prepare an initial mesh start¬ 
ing from a very coarse spatial mesh and a uniform partition of the time interval 
[0,1] with the step-size T = 0.05 until the user defined tolerances ttol = 10“^ and 
stol+ = 3 X 10“"^ are satisfied. Fig. T6|[T7 show the adaptive meshes and concen¬ 
trations at / = 0.3 and ^ = 1, where we can clearly see the flow-focusing due to the 
high-permeability. 



Figure 15: Example [7.4[ Permeability field (left) and velocity streamlines (right) 


Time vs DoFs and time vs time step-size plots are given in Fig.[^ We see that 
initially small time steps are used and then it reaches a steady time step, Fig. [T^ on 
the right. The number of DoFs increases (refinement dominates coarsening) mono- 
tonically after the meet of first high-permeability rock until the meet of third high- 
permeability rock and then the increase stops. Fig. on the left. This is mean¬ 
ingful since there is no sharp flow canalization after the third high-permeability 
rock. 
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t=0.3 


t=0.3 



Figure 16: Example [7.4[ Adaptive mesh (left) and concentration (right) at ^ = 0.3 



8 Conclusion 

We implemented a time-space adaptive algorithm for non-linear ADR equations 
based on utilizing the elliptic reconstruction technique to be able to use the elliptic 
a posteriori error estimator for the convection dominated parabolic problems with 
non-linear reaction mechanisms. We derived a posteriori error estimator in the 
Lr{L?) norm using backward Euler in time and SIPG in space. We 

demonstrated the performance of the algorithm by numerical studies. 
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Figure 18: Example [ t^ Evolution of DoFs (left) and time-steps T (right) 
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